import SimpleITK as sitk
import numpy as np
import matplotlib.pyplot as plt

# 读取 T1 图像和脑区分割图
t1_image = sitk.ReadImage('C:/Users/Administrator/Desktop/MNI152_T1_1mm.nii.gz')  # T1 加权 MRI 图像
seg_image = sitk.ReadImage('C:/Users/Administrator/Desktop/ID03_MOTOR_ICA.nii.gz')  # 脑区分割文件
print(t1_image)
# 获取 T1 图像的空间分辨率
t1_spacing = t1_image.GetSpacing()

# 重采样分割图像，使其和 T1 图像大小一致
resampler = sitk.ResampleImageFilter()
resampler.SetSize(t1_image.GetSize())
resampler.SetOutputSpacing(t1_spacing)
resampler.SetOutputOrigin(t1_image.GetOrigin())
resampler.SetOutputDirection(t1_image.GetDirection())
resampler.SetDefaultPixelValue(0)  # 如果在重采样时出现空值，使用 0 填充

# 使用线性插值进行重采样
resampler.SetInterpolator(sitk.sitkNearestNeighbor)  # 分割图通常用邻近插值
resampled_seg_image = resampler.Execute(seg_image)

# 将图像转换为 numpy 数组
t1_array1 = sitk.GetArrayFromImage(t1_image)
seg_array1 = sitk.GetArrayFromImage(resampled_seg_image)

t1_array = np.flip(t1_array1, axis=0)  # 例如，翻转 X 轴
seg_array = np.flip(seg_array1, axis=0)  # 例如，翻转 X 轴
# 选择一个切片进行可视化
slice_index = t1_array.shape[0] // 2  # 选择中间的切片

plt.subplot(2, 2, 1)
plt.imshow(t1_array[slice_index,:,:], cmap='gray')
plt.imshow(seg_array[slice_index,:,:], alpha=0.5, cmap='jet',vmin=-1.5, vmax=2)  # 使用透明度叠加分割图
plt.title('T1 Image with Segmentation')
plt.axis('off')


plt.subplot(2, 2, 2)
plt.imshow(t1_array[:,slice_index,:], cmap='gray')
plt.imshow(seg_array[:,slice_index,:], alpha=0.5, cmap='jet',vmin=-1.5, vmax=2)  # 使用透明度叠加分割图
plt.title('T1 Image with Segmentation')
plt.axis('off')

plt.subplot(2, 2, 3)
plt.imshow(t1_array[:,:,slice_index], cmap='gray')
plt.imshow(seg_array[:,:,slice_index], alpha=0.5, cmap='jet',vmin=-1.5, vmax=2)  # 使用透明度叠加分割图
plt.title('T1 Image with Segmentation')
plt.axis('off')

plt.show()